
`ifndef FLOAT_ARITH_PKG_SVH
 `define FLOAT_ARITH_PKG_SVH


package float_arith;

   
   typedef int unsigned uint32_t;
   typedef longint unsigned uint64_t;


   typedef enum    logic [1:0] {
				RN = 2'b00,
				RD = 2'b01,
				RU = 2'b10,
				RZ = 2'b11 
				} rounding_mode_t;
   
   
typedef struct packed {
   // logic       flush_to_zero;
   rounding_mode_t rounding_mode;
   // logic       overflow_en;	-- no useful ??
   // logic       underflow_en;
   logic       tininess_detect;	// 0:after rounding; 1:before rounding
} fp_cfg_t;

typedef struct packed {
   logic       sign;
   logic [7:0] exp;
   logic [22:0] sig;
} float32_t;

typedef struct packed {
   logic       sign;
   logic [10:0] exp;
   logic [51:0] sig;
} float64_t;


typedef struct packed {
   logic       invalid;
   logic       divbyzero;
   logic       overflow;
   logic       underflow;
   logic       inexact;
} excp_t;
   


typedef struct packed {
   bit       sign;
   bit [7:0] exp;
   bit [22:0] sig;
} float32_bit_t;

   
// typedef struct packed {
//    float32_t data;
//    logic [4:0] excp;
// } float32_result_t;


// typedef struct packed {
//    logic       hidden;
//    logic [22:0] sig;
// } float32_frac_t;
   
// typedef struct packed {
//    logic       carry;
//    // logic [22+1:0] frac;		// sig + hidden bit
//    float32_frac_t frac;
//    logic 	roundbit;
//    logic 	stickybit;
// } float32_mantissa_t;

typedef struct packed {
   // logic       carry;
   logic [22+1:0] frac;		// hidden bit + sig

   logic 	guard;
   logic 	round;
   logic 	sticky;
   // logic 	roundbit;
   // logic 	roundext;	// only for float32_sub
   // logic 	stickybit;
} float32_mantissa_t;

typedef struct packed {
   // logic       carry;
   logic [51+1:0] frac;		// hidden bit + sig

   logic 	guard;
   logic 	round;
   logic 	sticky;
   // logic 	roundbit;
   // logic 	roundext;	// only for float32_sub
   // logic 	stickybit;
} float64_mantissa_t;

   
// typedef struct packed {
//    logic       sign;
//    logic [7:0] exp;
//    float32_mantissa_t manti;   
// } float32_inter_t;		// intermediate type

typedef struct packed {
   logic       sign;
   logic [7:0] exp;		// unbiased exp
   float32_mantissa_t mantissa;
   
   // mantissa
   // logic       hidden;
   // logic [22:0] sig;
   // logic [22+1:0] frac;
   // logic 	roundbit;
   // logic 	stickybit;
} float32_inter_t;		// intermediate type

   

// typedef struct packed {
//    logic       carry;
//    logic       hiddenbit;
//    logic [22:0] sig;
//    logic 	roundbit;
//    logic 	stickybit;
// } float32_mantissa_t;

   
   const float32_t float32_default_nan = {1'b1, 8'hff, 23'h40_0000};
   const float64_t float64_default_nan = {1'b1, 11'h7ff, 52'h8_0000_0000_0000};

   
   
   function automatic float32_t pack_float32( input logic sign,
						     input logic [7:0] exp,
					      input logic [22:0] sig );

      return {sign, exp, sig};
   endfunction // pack_float32
   
					   
   
//    function automatic void float32_add( output float32_t result,
// 					output excp_t excp,
// 					input fp_cfg_t fp_cfg,
// 					input float32_t a, b );
//       //float32_t sum;

//       float32_t fa, fb;		// fa is the bigger one from a and b
      
//       logic [7+1:0] 					  exp_diff;
      
//       logic [22+1:0] 					  a_frac, b_frac;
//       logic [22+1:0] 					  tmp_frac;

      
//       logic [1:0] 					  round_bits; // for rounding operation
      
//       // float32_inter_t fa, fb;      
//       // float32_inter_t sum;
//       float32_t sum;

//       float32_mantissa_t ma, mb, msum;

//       excp = '0;

      
      
//       // 1. pre-process

//       // 2. shift  -> get sum.exp
//       // if(a.exp == 0)
//       // 	exp_diff = a.exp;
//       // else
//       // 	exp_diff = ;

//       sum.sign = a.sign;	// the sign of a & b must be the same 
      
//       exp_diff = a.exp - b.exp;
//       if($signed(exp_diff)>=0)		// a.exp > b.exp, b need shift right for exp_diff bits
// 	begin
// 	   fa = a;
// 	   fb = b;
// 	end
//       else
// 	begin
// 	   fa = b;
// 	   fb = a;
// 	   exp_diff = -$signed(exp_diff);
// 	end // else: !if($signed(exp_diff)>=0)

// 	  if(fa.exp && (fb.exp == 0)) // the smaller one is a subnormal
// 		exp_diff -= 1;				// the relative exponent is 1 not 0
	  

// 	  sum.exp = fa.exp;

//       ma = '{ carry:1'b0,
// 	  		  frac: {(fa.exp ? 1'b1 : 1'b0), fa.sig},
// 	  		  roundbit: 1'b0,
// 	  		  stickybit: 1'b0 };

//       mb = '{ carry:1'b0,
// 	  		  frac: {(fb.exp ? 1'b1 : 1'b0), fb.sig},
// 	  		  roundbit: 1'b0,
// 	  		  stickybit: 1'b0 };


//       // round bit
//       if(exp_diff && (exp_diff <= 24)) // 1<= exp_diff <= 24		--- FIXME: roundbit could be together with the shift of frac ???
// 		mb.roundbit = mb.frac[exp_diff-1];
//       else
// 		mb.roundbit = 0;

//       // calc sticky bit
// 	  begin
// 		 logic [23:0] 						  stickybit_mask;
// 		 // if(exp_diff == 0)
// 		 if(exp_diff <= 1)
// 		   stickybit_mask = '0;
// 		 else if(exp_diff <= 24)
// 		   stickybit_mask = (1<<(exp_diff-1))-1;
// 		 else
// 		   stickybit_mask = {24{1'b1}};

// 		 mb.stickybit = |(mb.frac & stickybit_mask);
// 	  end
	  
//       // if(exp_diff <= 23)
//       //   mb.frac >>= exp_diff;
//       // else
//       //   mb.frac = '0;
//       mb.frac >>= exp_diff;		// {mb.frac, mb.roundbit} = {1'b0, mb.frac} >> exp_diff;
	   
      
//       // 3. calc sum
//       msum = ma + mb;
      

//       // 4. normalize
//       if(msum.carry)
// 	begin
// 	   msum = {1'b0, msum.carry, msum.frac, msum.roundbit|msum.stickybit}; // shift right for one bit

// 	   // sum.mantissa = fsum;
// 	   sum.exp += 1;	// now sum.exp may be 8'hff
// 	end

//       // subnormal + subnormal and overflow to a normal
//       if( (fa.exp == 0) && (fb.exp == 0) && msum.frac[23] )
// 	sum.exp += 1;

      
//       // 5. check inexact/overflow/underflow
//       excp.inexact = msum.roundbit | msum.stickybit;

      
//       // if(sum.exp == 8'hff)
//       // 	begin
//       // 	   $display("sum.exp = 8'hff @", $time);
//       // 	end

      
// // ---------- round and pack ----------

      
//       // 6. round	--- but for subnormal num, first normalize or not ????
//       begin
// 	 // ==================== method 1 ====================	 
// 	 // logic [1:0] round_bits;
// 	 logic 	     round_nearest_even;

// 	 round_nearest_even = 0;
// 	 unique case(fp_cfg.rounding_mode)
// 	   RZ: round_bits = 2'b00;
// 	   RN: begin
// 	      round_bits = 2'b10;
// 	      round_nearest_even = {msum.roundbit, msum.stickybit} == 2'b10;
// 	   end
// 	   RU: round_bits = sum.sign ? 2'b00 : 2'b11;
// 	   RD: round_bits = sum.sign ? 2'b11 : 2'b00;
// 	 endcase // unique case (fp_cfg.rounding_mode)

// 	 msum += round_bits;
// 	 msum.frac[0] &= ~round_nearest_even;
	 

// 	 // // ==================== method 2 ====================	 
// 	 // logic [2:0] round_bits;
// 	 // unique case(fp_cfg.rounding_mode)
// 	 //   RZ: round_bits = 3'b000;
// 	 //   RN: begin
// 	 //      if({msum.roundbit, msum.stickybit}==2'b10) // in the case of a tie
// 	 // 	round_bits = {msum.frac[0], 2'b00};
// 	 //      else
// 	 // 	round_bits = 3'b010;
// 	 //   end
// 	 //   RU: round_bits = sum.sign ? 3'b000 : 3'b011;
// 	 //   RD: round_bits = sum.sign ? 3'b011 : 3'b000;
// 	 // endcase // unique case (fp_cfg.rounding_mode)

// 	 // msum += round_bits;
//       end
      
//       // normalize again	--- (with step 4) the two msum.carry occur once for max
//       if(msum.carry)
// 	begin
// 	   msum = {1'b0, msum.carry, msum.frac, msum.roundbit|msum.stickybit}; // shift right for one bit (the low 2 bits is not needed !!!)
// 	   sum.exp += 1;
// 	end
      
//       // repeat 4 and 5 ???
      
//       // check overflow
//       if(sum.exp == 8'hff)
// 	begin
// 	   excp.overflow = 1'b1;
// 	   excp.inexact = 1'b1;
	   
// 	   if(round_bits)	// the result is +-INF
// 	     begin
// 		sum.sig = '0;
// 	     end
// 	   else			// the result is +-max normal
// 	     begin
// 		sum.exp = 8'hfe;
// 		sum.sig = '1;
// 	     end
// 	end
//       else
// 	sum.sig = msum.frac[22:0];

//       result = sum;
//    endfunction // float32_add
   

   function automatic int count_leading_zeros(float32_mantissa_t mantissa);
      for(int i=$bits(mantissa)-1;i>=0;i--)
	if(mantissa[i] == 1'b1)
	  return ($bits(mantissa)-1-i);
      return $bits(mantissa);
   endfunction // count_leading_zeros
   
	  

   
   function automatic void float32_addsub( 
					   // for debug
					   output float32_mantissa_t o_mbig,
					   output float32_mantissa_t o_msml,
					   output bit [1:0] o_round_bits,
					   
					   output float32_t result,
					   output excp_t excp,
					   input fp_cfg_t fp_cfg,
					   input float32_t a, b,
					   input bit is_sub // 0:add, 1:sub
					   );

      float32_t fa, fb;		// fa is the bigger one from a and b
      
      logic [7+1:0] 					  exp_diff;
      logic [1:0] 					  round_bits; // for rounding operation
      logic 						  carry;
      
      int 						  cnt_zeros;
      logic 						  op_diff_sign;
      logic 						  op_all_subn;
      
      float32_t fbig, fsml;
      float32_mantissa_t mbig, msml;
      float32_inter_t sum;

      logic 						  a_is_qNaN, b_is_qNaN;
      logic 						  a_is_sNaN, b_is_sNaN;
      logic 						  a_is_NaN, b_is_NaN;
      logic 						  a_is_INF, b_is_INF;
      
      
      excp = '0;

      // =============== pre-process ===============
      a_is_qNaN = (a.exp==8'hff) && a.sig[22]; // the default qNaN ???
      b_is_qNaN = (b.exp==8'hff) && b.sig[22]; // or (&b.exp) && b.sig[22];

      a_is_sNaN = (a.exp==8'hff) && ~a.sig[22] && a.sig[21:0];
      b_is_sNaN = (b.exp==8'hff) && ~b.sig[22] && b.sig[21:0];

      a_is_NaN = (a.exp==8'hff) && a.sig;
      b_is_NaN = (b.exp==8'hff) && b.sig;

      a_is_INF = (a.exp==8'hff) && a.sig=='0;
      b_is_INF = (b.exp==8'hff) && b.sig=='0;


      // -------------------- process the NaN --------------------
      if(a_is_sNaN || b_is_sNaN)
	begin
	   excp.invalid = 1'b1;
	end
      
      if(a_is_sNaN)
	a.sig[22] = 1'b1;	// quiet a 

      if(b_is_sNaN)
	b.sig[22] = 1'b1;	// quiet b

      // a is sNaN
      if(a_is_sNaN && b_is_qNaN)
	begin
	   result = b;
	   return;
	end

      if(a_is_sNaN && ~b_is_NaN)
	begin
	   result = a;
	   return;
	end

      // a is qNaN
      if(a_is_qNaN && ~b_is_qNaN)
	begin
	   result = a;
	   return;
	end
      
      // a is not NaN
      if(~a_is_NaN && b_is_NaN)
	begin
	   result = b;
	   return;
	end
      
      if( (a_is_sNaN && b_is_sNaN) || (a_is_qNaN && b_is_qNaN))
	begin
	   // return the one which significand is larger
	   if(a.sig > b.sig)
	     result = a;
	   else if(a.sig < b.sig)
	     result = b;
	   else
	     result = a.sign ? b : a; // return the positive one
	   
	   return;
	end
      
      


      

      b.sign = b.sign ^ is_sub;	// convert a-b to a+(-b)
      op_diff_sign = a.sign ^ b.sign;
      
      // -------------------- process the INF -------------------- 
      if(a_is_INF && b_is_INF && op_diff_sign)
	begin
	   excp.invalid = 1'b1;
	   result = float32_default_nan;
	   return;
	end

      if(a_is_INF)
	begin
	   result = a;
	   return;
	end

      if(b_is_INF)
	begin
	   result = b;
	   return;
	end
      

      
      // op_diff_sign = a.sign ^ b.sign;
      
      
      // =============== compare exponents (irrelative with the sign of each op) =============== 
      exp_diff = a.exp - b.exp;
      if($signed(exp_diff)>=0)		// a.exp > b.exp, b need shift right for exp_diff bits
	begin
	   fbig = a;
	   fsml = b;
	end
      else
	begin
	   fbig = b;
	   fsml = a;
	   exp_diff = -$signed(exp_diff);
	end // else: !if($signed(exp_diff)>=0)

      // if(fbig.exp == 8'hff)
      // 	begin
      // 	   if(fbig.sig)		// the big one is NaN
      // 	     return xxx;
      // 	   else
      // 	     begin
      // 		if(op_diff_sign) // -INF+(+INF) or +inf + (-inf) => qNaN, invalid
      // 		  begin
      // 		     excp.invalid = 1'b1;
      // 		     return qNaN;
      // 		  end
      // 		else
      // 		  return fbig;	// the big one is INF, but the sign ???
      // 	     end
      // 	end

      // if()

      

      
      if(fbig.exp && (fsml.exp == 0)) // the smaller one is a subnormal
	exp_diff -= 1;				// the relative exponent is 1 not 0	--- FIXME: or first shift left the subnormal num for one bit ???
	  
      sum.exp = fbig.exp;

 
      // =============== mantissa of fbig/fsml ===============
      mbig = '{ // carry: 1'b0,
	       frac: {(fbig.exp ? 1'b1 : 1'b0), fbig.sig},
	       guard:1'b0,
	       round:1'b0,
	       sticky:1'b0 };

      msml.frac = {(fsml.exp ? 1'b1 : 1'b0), fsml.sig};
      // calc sticky bit
      begin
	 logic [23:0] 						  sticky_mask;
	 if(exp_diff <= 2)
	   sticky_mask = '0;
	 else if(exp_diff <= 25) // 3 ~ 25
	   sticky_mask = (1<<(exp_diff-2))-1;
	 else
	   sticky_mask = {24{1'b1}};
	 
	 msml.sticky = |(msml.frac & sticky_mask);
      end
      // cal guard, round bit
      {msml.frac, msml.guard, msml.round} = {msml.frac, 2'b00} >> exp_diff;
	   

      o_mbig = mbig;
      o_msml = msml;
      
      // =============== calc intermediate result (the sign is dealed) ===============
      op_all_subn = {mbig.frac[23], msml.frac[23]} == 2'b00;
      if(op_diff_sign)
	begin
      	   if(fbig.sign)
      	     mbig = ~mbig;
	   else
	     msml = ~msml;
	end
      
      {carry, sum.mantissa} = mbig + msml + op_diff_sign; // if the sum/diff is 0, then carry must be 1, 
      
      sum.sign = op_diff_sign ? ~carry : fbig.sign; // if the diff is 0, the sign is 0 --- but need according to the round mode ??
      if(op_diff_sign && ~carry)	// the result is negative
	sum.mantissa = ~sum.mantissa + 1'b1;

      // normalize(or un-normalize for subnormal) -> round -> normalize -> pack
      if(~op_diff_sign)		 // addition
	begin
	   if(carry)		 // frac overflow, shift right for 1 bit
	     sum.mantissa = {1'b1, sum.mantissa.frac, sum.mantissa.guard, sum.mantissa.round | sum.mantissa.sticky};
	   if(carry || (op_all_subn && sum.mantissa.frac[23])) // 0.xxx + 0.xxx => 1.xxx
	     // sum_exp_increase = 1'b1;
	     sum.exp++;
	end
      
      if(sum.mantissa.frac[23]==0) // need shift left	--- only subtraction ???
	begin
	   // count leading zeros
	   cnt_zeros = count_leading_zeros(sum.mantissa);

	   if(sum.exp > cnt_zeros)
	     begin
		sum.mantissa <<= cnt_zeros; // the sticky bit is ??
		sum.exp -= cnt_zeros;
	     end
	   else if(sum.exp)	// 1 <= sum.exp <= cnt_zeros
	     begin
		sum.mantissa <<= sum.exp-1;
		sum.exp = 0;
	     end
	   // else			// sum.exp == 0, not shift

	   // ==========
	   if(sum.mantissa.frac == 0) // the result is exactly 0
	     begin
		sum.exp = 0;
		sum.sign = (fp_cfg.rounding_mode==RD); // the result depends on the rounding mode when a,b have different sign
	     end
	   // ==========	
	   
	end // if (sum.mantissa.frac[23]==0)
      

      // guard bit is the actual round bit
      sum.mantissa.sticky = sum.mantissa.round | sum.mantissa.sticky;
      
      // 5. check inexact/overflow/underflow
      excp.inexact = sum.mantissa.guard | sum.mantissa.sticky;

      
      // if(sum.exp == 8'hff)
      // 	begin
      // 	   $display("sum.exp = 8'hff @", $time);
      // 	end
      
// ---------- round and pack (the input is sum) ----------
      
      // 6. round	--- but for subnormal num, first normalize or not ????
      begin
	 // ==================== method 1 ====================	 
	 // logic [1:0] round_bits;
	 logic 	     round_nearest_even;

	 round_nearest_even = 0;
	 unique case(fp_cfg.rounding_mode)
	   RZ: round_bits = 2'b00;
	   RN: begin
	      round_bits = 2'b10;
	      round_nearest_even = {sum.mantissa.guard, sum.mantissa.sticky} == 2'b10;
	   end
	   RU: round_bits = sum.sign ? 2'b00 : 2'b11;
	   RD: round_bits = sum.sign ? 2'b11 : 2'b00;
	 endcase // unique case (fp_cfg.rounding_mode)


	 // for debug
	 o_round_bits =  {sum.mantissa.guard, sum.mantissa.sticky};
	 
	 
	 carry = 0;
	 {carry, sum.mantissa.frac, sum.mantissa.guard, sum.mantissa.sticky} += round_bits;
										     
	 sum.mantissa.frac[0] &= ~round_nearest_even;

	 
	 // // ==================== method 2 ====================	 
	 // logic [2:0] round_bits;
	 // unique case(fp_cfg.rounding_mode)
	 //   RZ: round_bits = 3'b000;
	 //   RN: begin
	 //      if({msum.round, msum.sticky}==2'b10) // in the case of a tie
	 // 	round_bits = {msum.frac[0], 2'b00};
	 //      else
	 // 	round_bits = 3'b010;
	 //   end
	 //   RU: round_bits = sum.sign ? 3'b000 : 3'b011;
	 //   RD: round_bits = sum.sign ? 3'b011 : 3'b000;
	 // endcase // unique case (fp_cfg.rounding_mode)

	 // msum += round_bits;

	 // // ==================== method 3 ====================
	 //    logic 	     increment;
	 //    unique case(fp_cfg.rounding_mode)
	 //      RZ: begin
	 // 	 increment = 0;
	 //      end
	 //      RN: begin
	 // 	 unique case({sum.mantissa.round, sum.mantissa.sticky})
	 // 	   2'b00:
	 // 	     increment = 0;
	 // 	   2'b01:
	 // 	     increment = 0;
	 // 	   2'b10:		// in the case of a tie
	 // 	     increment = sum.frac[0];
	 // 	   2'b11:
	 // 	     increment = 1'b1;
	 // 	 endcase // unique case ({sum.round, sum.sticky})
	 //      end
	 //      RU: increment = sum.sign ? 1'b0 : ({1'b0, sum.mantissa.round, sum.mantissa.sticky} + 3'b011)[2];
	 //      RD: increment = sum.sign ? ({1'b0, sum.mantissa.round, sum.mantissa.sticky} + 3'b011)[2] : 1'b0;
	 //    endcase // unique case (fp_cfg.rounding_mode)

	 //    {carry, sum.mantissa.frac} += increment;
      end
      
      // normalize again	--- (with step 4) the two msum.carry occur once for max
      if(carry)
	begin
	   sum.mantissa >>= 1;
	   // sum = {1'b0, msum.carry, msum.frac, msum.round|msum.sticky}; // shift right for one bit (the low 2 bits is not needed !!!)
	   sum.exp += 1;
	end
      


      // repeat 4 and 5 ???
      
      // check overflow
      if(sum.exp == 8'hff)
      	begin
      	   excp.overflow = 1'b1;
      	   excp.inexact = 1'b1;
	   
      	   if(round_bits)	// the result is +-INF
      	     begin
		result.exp = 8'hff;
      		result.sig = '0;
      	     end
      	   else			// the result is +-max normal
      	     begin
      		result.exp = 8'hfe;
      		result.sig = '1;
      	     end
      	end
      else
	begin
	   result.exp = sum.exp;
      	   result.sig = sum.mantissa.frac[22:0];
	end

      result.sign = sum.sign;
   endfunction // float32_sub
   
   

//    function automatic void float32_sub_test( output float32_t result,
// 					output excp_t excp,
// 					input fp_cfg_t fp_cfg,
// 					input float32_t a, b );
//       //float32_t sum;

//       // float32_t fa, fb;		// fa is the bigger one from a and b
      
//       logic [7+1:0] 					  exp_diff;
      
//       logic [22+1:0] 					  a_frac, b_frac;
//       logic [22+1:0] 					  tmp_frac;

      
//       logic [1:0] 					  round_bits; // for rounding operation

//       float32_t fbig, fsml;
      
      
//       float32_inter_t fa, fb;  
//       float32_inter_t sum;
//       // float32_t sum;
      

//       float32_mantissa_t ma, mb, msum;

//       excp = '0;


//       // a + b  // a is positive, b is negative, but the sign of result ???
      
      
//       // 1. pre-process
//       // sum.sign = a.sign;	// the sign of a & b must be the same 
      
//       exp_diff = a.exp - b.exp;
//       if($signed(exp_diff)>=0)		// a.exp > b.exp, b need shift right for exp_diff bits
// 	begin
// 	   fbig = a;
// 	   fsml = b;
// 	end
//       else
// 	begin
// 	   fbig = b;
// 	   fsml = a;
// 	   exp_diff = -$signed(exp_diff);
// 	end // else: !if($signed(exp_diff)>=0)

//       if(fbig.exp && (fsml.exp == 0)) // the smaller one is a subnormal
// 	exp_diff -= 1;				// the relative exponent is 1 not 0
	  

      
//       sum.exp = fbig.exp;

//       fa = '{ sign:fbig.sign,
//       	      exp:fbig.exp,
// 	      frac: {(fbig.exp ? 1'b1 : 1'b0), fbig.sig},
//       	      roundbit: 1'b0,
//       	      stickybit: 1'b0 }; 

//       fb = '{ sign:fsml.sign,
//       	      exp:fsml.exp,
// 	      frac: {(fsml.exp ? 1'b1 : 1'b0), fsml.sig},
//       	      roundbit: 1'b0,
//       	      stickybit: 1'b0 };

      
//       // ma = '{ carry:1'b0,
//       // 	      frac: {(fbig.exp ? 1'b1 : 1'b0), fbig.sig},
//       // 	      roundbit: 1'b0,
//       // 	      stickybit: 1'b0 };
      
//       // mb = '{ carry:1'b0,
//       // 	      frac: {(fsml.exp ? 1'b1 : 1'b0), fsml.sig},
//       // 	      roundbit: 1'b0,
//       // 	      stickybit: 1'b0 };


//       // round bit
//       if(exp_diff && (exp_diff <= 24)) // 1<= exp_diff <= 24
// 	fb.roundbit = fb.frac[exp_diff-1];
//       else
// 	fb.roundbit = 0;

//       // calc sticky bit
//       begin
// 	 logic [23:0] 						  stickybit_mask;
// 	 if(exp_diff == 0)
// 	   stickybit_mask = '0;
// 	 else if(exp_diff <= 24)
// 	   stickybit_mask = (1<<(exp_diff-1))-1;
// 	 else
// 	   stickybit_mask = {24{1'b1}};
	 
// 	 fb.stickybit = |(fb.frac & stickybit_mask);
//       end
	  
//       // if(exp_diff <= 23)
//       //   mb.frac >>= exp_diff;
//       // else
//       //   mb.frac = '0;
//       fb.frac >>= exp_diff;
	   
      
//       // 3. calc sum
//       // msum = ma + mb;
//       begin
// 	 logic carry;
//       {carry, sum.frac, sum.roundbit, sum.stickybit} = {1'b0, fa.frac, fa.roundbit, fa.stickybit} -
// 						       {1'b0, fb.frac, fb.roundbit, fb.stickybit};
// end
//       // ==================================================
//       //{msum.frac, msum.roundbit, msum.stickybit} = ;
//       sum.sign = msum.carry;
//       if(msum.carry)		// the borrow bit
// 	begin
// 	   msum = -msum;	// the magnitude of msum
// 	end

      
//       // normalize(or un-normalize for subnormal) -> round -> normalize -> pack
//       // ==================================================      

      
//       // 4. normalize
//       if(msum.carry)		// the borrow bit
// 	begin
// 	   msum = {1'b0, msum.carry, msum.frac, msum.roundbit|msum.stickybit}; // shift right for one bit

// 	   // sum.mantissa = fsum;
// 	   sum.exp += 1;	// now sum.exp may be 8'hff
// 	end

//       // subnormal + subnormal and overflow to a normal
//       if( (fa.exp == 0) && (fb.exp == 0) && msum.frac[23] )
// 	sum.exp += 1;

      
//       // 5. check inexact/overflow/underflow
//       excp.inexact = msum.roundbit | msum.stickybit;

      
//       // if(sum.exp == 8'hff)
//       // 	begin
//       // 	   $display("sum.exp = 8'hff @", $time);
//       // 	end

      
// // ---------- round and pack ----------

      
//       // 6. round	--- but for subnormal num, first normalize or not ????
//       begin
// 	 // ==================== method 1 ====================	 
// 	 // logic [1:0] round_bits;
// 	 logic 	     round_nearest_even;

// 	 round_nearest_even = 0;
// 	 unique case(fp_cfg.rounding_mode)
// 	   RZ: round_bits = 2'b00;
// 	   RN: begin
// 	      round_bits = 2'b10;
// 	      round_nearest_even = {msum.roundbit, msum.stickybit} == 2'b10;
// 	   end
// 	   RU: round_bits = sum.sign ? 2'b00 : 2'b11;
// 	   RD: round_bits = sum.sign ? 2'b11 : 2'b00;
// 	 endcase // unique case (fp_cfg.rounding_mode)

// 	 msum += round_bits;
// 	 msum.frac[0] &= ~round_nearest_even;
	 

// 	 // // ==================== method 2 ====================	 
// 	 // logic [2:0] round_bits;
// 	 // unique case(fp_cfg.rounding_mode)
// 	 //   RZ: round_bits = 3'b000;
// 	 //   RN: begin
// 	 //      if({msum.roundbit, msum.stickybit}==2'b10) // in the case of a tie
// 	 // 	round_bits = {msum.frac[0], 2'b00};
// 	 //      else
// 	 // 	round_bits = 3'b010;
// 	 //   end
// 	 //   RU: round_bits = sum.sign ? 3'b000 : 3'b011;
// 	 //   RD: round_bits = sum.sign ? 3'b011 : 3'b000;
// 	 // endcase // unique case (fp_cfg.rounding_mode)

// 	 // msum += round_bits;
//       end
      
//       // normalize again	--- (with step 4) the two msum.carry occur once for max
//       if(msum.carry)
// 	begin
// 	   msum = {1'b0, msum.carry, msum.frac, msum.roundbit|msum.stickybit}; // shift right for one bit (the low 2 bits is not needed !!!)
// 	   sum.exp += 1;
// 	end
      
//       // repeat 4 and 5 ???
      
//       // check overflow
//       if(sum.exp == 8'hff)
// 	begin
// 	   excp.overflow = 1'b1;
// 	   excp.inexact = 1'b1;
	   
// 	   if(round_bits)	// the result is +-INF
// 	     begin
// 		sum.sig = '0;
// 	     end
// 	   else			// the result is +-max normal
// 	     begin
// 		sum.exp = 8'hfe;
// 		sum.sig = '1;
// 	     end
// 	end
//       else
// 	sum.sig = msum.frac[22:0];


      
//       result = sum;
//    endfunction // float32_add


   
   // function automatic int func_ref(input logic x1,x2, r, s, R1, R0);
   //    int 	  res;
   //    logic [4:0] sum1, sum2;
   //    sum1 = {x1, 2'b00} + {x2, r, s};
   //    sum2 = {sum1[2], |sum1[1:0], 1'b0} + {R1, R0, 1'b0};
   //    res = sum1[4:3] + sum2[4:3];
   //    return res;
   // endfunction // func_2


   // typedef logic [3:0] queue_bits_t[$];

   
   // function automatic queue_bits_t func_1(input logic x1,x2, r, s, input int correct_res);
   //    int res;
   //    queue_bits_t solutions;
      
   //    for(int i=0;i<=4'b1111;i++)
   // 	begin
   // 	   res = {x1, r, s} + {x2, 2'b00} + i[3:0];
   // 	   res >>= 3;
   // 	   if(res === correct_res)
   // 	     solutions.push_back(i[3:0]);
   // 	end
   //    return solutions;
   // endfunction // func_1
   
   
   // function automatic void test();
   //    int 	  res;
   //    int 	  ia, ib;
   //    logic 	  x1,x2, r, s, R2, R1, R0;
   //    for(int i=0;i<=6'b111111;i++)
   // 	begin
   // 	   queue_bits_t solutions;
	   
   // 	   {x1, x2, r, s, R1, R0} = i[5:0];
   // 	   res = func_ref(x1,x2, r, s, R1, R0); // the correct val
	   
   // 	   solutions = func_1(x1, x2, r, s, res);
   // 	   $display("input:%b, solutions: ", i[5:0], solutions);
	   
   // 	   // ia = func_1(x1,x2, r, s, 1'b0, R1, R0);
   // 	   // ib = func_1(x1,x2, r, s, 1'b1, R1, R0);
   // 	   // $display("input=%b, result:%b, func_1:%b/%b, %s equal", i[5:0], res[2:0], ia[2:0], ib[2:0], ia===res ? "a" : ib===res ? "b" : "not");

	     
   // 	end
   // endfunction // test
   


   // floating-point multiply
   function automatic void float32_mul( 
					output float32_t result,
					output excp_t excp,
					input fp_cfg_t fp_cfg,
					input float32_t a, b, c
					);
      bit [9:0] 			      ex;
      bit [47:0] 			      prod, prod_sft, prod_t;
      bit [23:0] 			      frac_a, frac_b;
      bit 				      sign;
      // for rounding
      bit 				      round, sticky;
      bit [1:0] 			      round_inc_t, round_inc;
      bit 				      round_nearest_even;
      bit 				      round_carry;
      bit [1:0] 			      tmp_bits;
      // after rounding
      bit [23:0] 			      frac_prod_rnd;
      bit 				      frac_rnd_carry;
      // result
      bit [7:0] 			      result_exp;
      bit [22:0] 			      result_sig;


      sign = a.sign ^ b.sign;	// neight a nor b is 0
      
      frac_a = a.exp ? {1'b1, a.sig} : {a.sig, 1'b0};
      frac_b = b.exp ? {1'b1, b.sig} : {b.sig, 1'b0};
      prod = frac_a * frac_b;

      ex = a.exp + b.exp - 127 - 1;

      if(ex[9])			// ex < 0
	begin
// ============================== ex < 0 ==============================	   
	   prod_sft = prod >> (~ex[6:0]);
	   // round/sticky bits
	   if( (~ex[6:0]) > $bits(prod))
	     begin
		round = 1'b0;
		sticky = |prod;
	     end
	   else
	     begin
		bit [47:0] prod_t;
		prod_t = prod << ($bits(prod) - (~ex[6:0]));
		round = prod_t[47];
		sticky = |prod_t[46:0];
	     end // else: !if( (~ex[6:0]) > $bits(prod))

	   round_nearest_even = 0;
   	   unique case(fp_cfg.rounding_mode)
   	     RZ: round_inc_t = 2'b00;
   	     RN: begin
   		round_inc_t = 2'b10;
   		round_nearest_even = {round, sticky} == 2'b10;
   	     end
  	     RU: round_inc_t = {2{~sign}};
   	     RD: round_inc_t = {2{sign}};
   	   endcase // unique case (fp_cfg.rounding_mode)

	   round_inc = round_nearest_even ? {prod_sft[24], 1'b0} : round_inc_t;
	   {round_carry, tmp_bits} = {round, sticky} + round_inc;

	   {frac_rnd_carry, frac_prod_rnd} = prod_sft[47:24] + round_carry;

	   result_exp = {7'b0, frac_prod_rnd[23]};
	   result_sig = frac_prod_rnd[22:0];

// ============================== ex < 0 ==============================	   
	   
	   
	end // if (ex[9])
      else  // ex >= 0
	begin
	   bit [7:0] lz_n;
// ============================== ex >= 0 ==============================	   
	   // count leading zeros
	   lz_n = 0;
	   for(int i=47;i>=25;i--)
	     if(prod[i] == 1'b0)
	       lz_n++;
	     else
	       break;

	   if(ex[8:0] >= lz_n)
	     begin
		prod_sft = prod << lz_n;
		// 
		// result_exp = ex[8:0] + 2 - lz_n + frac_rnd_carry;
	     end
	   else
	     begin
		prod_sft = prod << (ex[7:0]+1);
		// result_exp = {6'b0, frac_rnd_carry, MSB};
	     end

	   round = prod_sft[23];
	   sticky = |prod_sft[22:0];
	   
	   round_nearest_even = 0;
   	   unique case(fp_cfg.rounding_mode)
   	     RZ: round_inc_t = 2'b00;
   	     RN: begin
   		round_inc_t = 2'b10;
   		round_nearest_even = {round, sticky} == 2'b10;
   	     end
  	     RU: round_inc_t = {2{~sign}};
   	     RD: round_inc_t = {2{sign}};
   	   endcase // unique case (fp_cfg.rounding_mode)

	   round_inc = round_nearest_even ? {prod_sft[24], 1'b0} : round_inc_t;
	   {round_carry, tmp_bits} = {round, sticky} + round_inc;

	   {frac_rnd_carry, frac_prod_rnd} = prod_sft[47:24] + round_carry;


	   result_exp = (ex[8:0] >= lz_n) ? ex+2-lz_n+frac_rnd_carry : {7'b0, frac_rnd_carry, prod_sft[47]};
	   
	   result_sig = frac_prod_rnd[22:0];

// ============================== ex >= 0 ==============================
	   
	end // else: !if(ex[9])
      

      result = {sign, result_exp, result_sig};
      excp = '0;

   endfunction
   
   
   // fused mul-add : o(a*b+c), round only once
   function automatic void float32_fma(
				       output float32_t result,
				       output excp_t excp,
				       input fp_cfg_t fp_cfg,
				       input float32_t a, b, c,
				       // flags for debug
				       output bit debug
				       );
      bit [9:0] 				  ex;
      bit [47:0] 				  prod, prod_sft, prod_t;
      bit [23:0] 				  frac_a, frac_b;
      bit 					  sign;
      bit [4:0] 				  lz_n;
      // for rounding
      bit 					  round, sticky;
      bit [1:0] 				  round_inc_t, round_inc;
      bit 					  round_nearest_even;
      bit 					  round_carry;
      bit [1:0] 				  tmp_bits;
      // after rounding
      bit [23:0] 				  frac_prod_rnd;
      bit 					  frac_rnd_carry;
      // result
      bit [7:0] 				  result_exp;
      bit [22:0] 				  result_sig;


      // variables for FMA
      bit [9:0] 				  diff_base;
      
      bit [23:0] 				  frac_c;
      bit [23:0] 				  frac_c_sft;
      bit [2:0] 				  frac_c_grs; // g/r/s bits
      bit 					  c_is_subn;
      bit [47:0] 				  frac_c_ext; // if c shift right
      
      bit [7:0] 				  expn; // intermediate exponent of the sum of FMA
      bit [23:0] 				  sum;
      bit [2:0] 				  sum_grs;
      bit [23:0] 				  sum_rnd;
      bit  					  sum_rnd_carry;
      bit 					  carry;

      bit [2:0] 				  prod_grs; // guard/round/sticky bits


      sign = a.sign ^ b.sign;	// neight a nor b is 0
      
      frac_a = a.exp ? {1'b1, a.sig} : {a.sig, 1'b0};
      frac_b = b.exp ? {1'b1, b.sig} : {b.sig, 1'b0};
      prod = frac_a * frac_b;

      ex = a.exp + b.exp - 127 - 1;

      // for FMA
      diff_base = c.exp - a.exp - b.exp + 126; //c.exp - ex - 2;
      c_is_subn = c.exp == '0;
      frac_c = { (c_is_subn ? 1'b0 : 1'b1), c.sig };


      // debug
      // debug = ex[9];
      // debug = ~diff_base[9];
      debug = 0;

      if(diff_base[9] == 1'b0)	// diff_base>=0
	begin
	   // prod_sft = prod >> (diff_base[8:0]+c_is_subn); // c is normal or subnormal
	   // expn = c_is_subn ? prod_sft[47] : c.exp; // this signal should be calculated quickly ??

	   // ==================== only for normal number ====================
	   prod_sft = prod >> diff_base[8:0];
	   frac_c_sft = frac_c;
	   expn = c.exp;

	   case(diff_base[8:0])
	     9'd0 : prod_grs = {prod[23], prod[22], |prod[21:0]};
	     9'd1 : prod_grs = {prod[24], prod[23], |prod[22:0]};
	     9'd2 : prod_grs = {prod[25], prod[24], |prod[23:0]};
	     9'd3 : prod_grs = {prod[26], prod[25], |prod[24:0]};
	     9'd4 : prod_grs = {prod[27], prod[26], |prod[25:0]};
	     9'd5 : prod_grs = {prod[28], prod[27], |prod[26:0]};
	     9'd6 : prod_grs = {prod[29], prod[28], |prod[27:0]};
	     9'd7 : prod_grs = {prod[30], prod[29], |prod[28:0]};
	     9'd8 : prod_grs = {prod[31], prod[30], |prod[29:0]};
	     9'd9 : prod_grs = {prod[32], prod[31], |prod[30:0]};
	     9'd10 : prod_grs = {prod[33], prod[32], |prod[31:0]};
	     9'd11 : prod_grs = {prod[34], prod[33], |prod[32:0]};
	     9'd12 : prod_grs = {prod[35], prod[34], |prod[33:0]};
	     9'd13 : prod_grs = {prod[36], prod[35], |prod[34:0]};
	     9'd14 : prod_grs = {prod[37], prod[36], |prod[35:0]};
	     9'd15 : prod_grs = {prod[38], prod[37], |prod[36:0]};
	     9'd16 : prod_grs = {prod[39], prod[38], |prod[37:0]};
	     9'd17 : prod_grs = {prod[40], prod[39], |prod[38:0]};
	     9'd18 : prod_grs = {prod[41], prod[40], |prod[39:0]};
	     9'd19 : prod_grs = {prod[42], prod[41], |prod[40:0]};
	     9'd20 : prod_grs = {prod[43], prod[42], |prod[41:0]};
	     9'd21 : prod_grs = {prod[44], prod[43], |prod[42:0]};
	     9'd22 : prod_grs = {prod[45], prod[44], |prod[43:0]};
	     9'd23 : prod_grs = {prod[46], prod[45], |prod[44:0]};
	     9'd24 : prod_grs = {prod[47], prod[46], |prod[45:0]};
	     9'd25 : prod_grs = {1'b0, prod[47], |prod[46:0]};
	     default:
	       prod_grs = {1'b0, 1'b0, |prod[47:0]};
	   endcase // case (diff_base[8:0])

	   frac_c_grs = 3'b000;

	   $finish;
	end
      else			// diff_base < 0
	begin
	   // count leading zeros
	   lz_n = 0;		
	   for(int i=47;i>=25;i--)
	     if(prod[i]==1'b0)
	       lz_n++;		// 0 ~ 23
	     else
	       break;
	   
	   // $display("lz_n = %d, @", lz_n, $time);
	   
	   if(~diff_base[8:0]>=lz_n)
	     begin
		bit [8:0] sft_n;
		bit [47+2:0] sft_ext;
		bit [47:0]   sum_ext;
		bit 	     sum_ext_tmp;
		bit [2:0]    sum_ext_grs;
		bit 	     sum_ext_carry;
		bit [47:0]   sticky_bits;
		
		sft_n = ~diff_base[8:0]+9'd1-lz_n;
		
		
	   	prod_sft = prod << lz_n;
		// prod_grs = {prod_sft[23], prod_sft[22], |prod_sft[21:0]};
	   	expn = ex + 2 - lz_n;

		frac_c_ext = {frac_c, 24'b0} >> sft_n;
		
		if(sft_n <= 48)
		  sticky_bits = {frac_c, 24'b0} << (48-sft_n);
		else
		  sticky_bits = |frac_c;
		  // sticky_bits = '0; // for debug	// also right ??? to simplify the logic ??

		{sum_ext_carry, sum_ext, sum_ext_tmp} = {prod_sft, 1'b0} + {frac_c_ext, |sticky_bits};
		
		if(sum_ext_carry)
		  begin
		     sum_ext_tmp = sum_ext_tmp | sum_ext[0];
		     sum_ext = sum_ext >> 1;
		     expn = expn + 1;
		  end
		
		round = sum_ext[23];
		sticky = |sum_ext[22:0] | sum_ext_tmp;
      
      round_nearest_even = 0;
      unique case(fp_cfg.rounding_mode)
   	RZ: round_inc_t = 2'b00;
   	RN: begin
   	   round_inc_t = 2'b10;
   	   round_nearest_even = {round, sticky} == 2'b10;
   	end
  	RU: round_inc_t = {2{~sign}};
   	RD: round_inc_t = {2{sign}};
      endcase // unique case (fp_cfg.rounding_mode)

      round_inc = round_nearest_even ? {sum_ext[24], 1'b0} : round_inc_t;
      {round_carry, tmp_bits} = {round, sticky} + round_inc;

      
      {sum_rnd_carry, sum_rnd} = sum_ext[47:24] + round_carry;

      
      result_sig = sum_rnd[22:0];
      if(sum_rnd_carry)
	result_exp = expn + 1;
      else
	result_exp = expn;
     
		
		
		
		result = {sign, result_exp, result_sig};
		excp = '0;


		return;
		
		//   9'd0 : frac_c_ext = ;
			    
		// // prod_grs = {prod_sft[23], prod_sft[22], |prod_sft[21:0]};
		// // frac_c_grs = 3'b000; // not the correct value ???
		
		// // xxx = frac_c << (24-sft_n);
		
		// case(sft_n)
		//   9'd0 : frac_c_sticky_bits = { 24'b0 };
		//   9'd1 : frac_c_sticky_bits = { 24'b0 };
		//   9'd2 : frac_c_sticky_bits = { 24'b0 };
		//   9'd3 : frac_c_sticky_bits = { 24'b0 };
		  
		//   9'd1 : frac_c_grs_x = { frac_c[0], 23'b0};
		//   9'd2 : frac_c_grs_x = { frac_c[1:0], 22'b0};
		//   9'd3 : frac_c_grs_x = { frac_c[2:0], 21'b0};
		//   9'd4 : frac_c_grs_x = { frac_c[3:0], 20'b0};
		
		//   default:;
		// endcase // case (sft_n)


		debug = 1;
		
	     end
	   else
	     begin
	   	prod_sft = prod << (~diff_base[8:0]+9'd1);
	   	expn = c.exp;
		prod_grs = {prod_sft[23], prod_sft[22], |prod_sft[21:0]};
		
		frac_c_sft = frac_c;	// not shift
		frac_c_grs = 3'b000;
	     end // else: !if(~diff_base[8:0]>=lz_n)
	end // else: !if(diff_base[9] == 1'b0)


      {carry, sum, sum_grs} = {prod_sft[47:24], prod_grs} + {frac_c_sft, frac_c_grs}; // first get the carry but not need the whole sum ??

      if(carry)
	begin
	   sum_grs = {sum[0], sum_grs[2], |sum_grs[1:0]};
	   sum = sum >> 1;
	   expn = expn + 1;
	end
      
      
      round = sum_grs[2];
      sticky = |sum_grs[1:0];
      
      round_nearest_even = 0;
      unique case(fp_cfg.rounding_mode)
   	RZ: round_inc_t = 2'b00;
   	RN: begin
   	   round_inc_t = 2'b10;
   	   round_nearest_even = {round, sticky} == 2'b10;
   	end
  	RU: round_inc_t = {2{~sign}};
   	RD: round_inc_t = {2{sign}};
      endcase // unique case (fp_cfg.rounding_mode)

      round_inc = round_nearest_even ? {sum[0], 1'b0} : round_inc_t;
      {round_carry, tmp_bits} = {round, sticky} + round_inc;

      
      {sum_rnd_carry, sum_rnd} = sum[23:0] + round_carry;

      
      result_sig = sum_rnd[22:0];
      if(sum_rnd_carry)
	result_exp = expn + 1;
      else
	result_exp = expn;
     
      result = {sign, result_exp, result_sig};
      excp = '0;

   endfunction // float32_fma





   // fused mul-add : o(a*b+c), round only once
   function automatic void float32_fma_new(
				       output float32_t result,
				       output excp_t excp,
				       input fp_cfg_t fp_cfg,
				       input float32_t a, b, c,
				       // flags for debug
				       output bit debug
				       );
      bit [9:0] 				  ex;
      bit [47:0] 				  prod, prod_sft, prod_t;
      bit [23:0] 				  frac_a, frac_b;
      bit 					  sign;
      bit [4:0] 				  lz_n;
      // for rounding
      bit 					  round, sticky;

      bit 					  round_carry;
      bit [1:0] 				  tmp_bits;
      // after rounding
      bit [23:0] 				  frac_prod_rnd;
      bit 					  frac_rnd_carry;
      // result
      bit [7:0] 				  result_exp;
      bit [22:0] 				  result_sig;


      // variables for FMA
      bit [9:0] 				  diff_base;
      
      bit [23:0] 				  frac_c;

      
      bit [2:0] 				  frac_c_grs; // g/r/s bits
      bit 					  c_is_subn;
      bit [47:0] 				  frac_c_ext; // if c shift right
      
      bit [7:0] 				  expn; // intermediate exponent of the sum of FMA
      bit [23:0] 				  sum;
      bit [2:0] 				  sum_grs;
      bit [23:0] 				  sum_rnd;
      bit  					  sum_rnd_carry;
      bit 					  carry;

      bit [2:0] 				  prod_grs; // guard/round/sticky bits


      // ==================== rounding (common for ??) ====================
      bit [3:0] 				  round_inc_orig, round_inc_t, round_inc;
      bit 					  lsb_bit;
      bit [1:0] 				  round_bits;
      bit 					  round_nearest_even;
      // logic       overflow_val_sel;



      // for branch 1
      	 bit [4:0] lzn;
	 
	 bit [47:0] frac_c_sft;
	 bit 	     frac_c_sticky;
	 
	 bit [9:0]   sft_bits;	// the range is ok ??
	 bit [47:0]  sum_diff_t, sum_diff_lzd;
	 bit 	     sum_diff_ext;

	 bit [27:0]  sum_diff, sum_diff_rnd;
	 
	 bit 	    is_true_sub;

      float32_t lzd_path_result;
      excp_t lzd_path_excp;

      bit 	    sel_lzd_path;


      sign = a.sign ^ b.sign;	// neight a nor b is 0
      
      frac_a = a.exp ? {1'b1, a.sig} : {a.sig, 1'b0};
      frac_b = b.exp ? {1'b1, b.sig} : {b.sig, 1'b0};
      prod = frac_a * frac_b;

      ex = a.exp + b.exp - 127 - 1;

      // for FMA
      diff_base = c.exp - a.exp - b.exp + 126; //c.exp - ex - 2;
      c_is_subn = c.exp == '0;
      frac_c = { (c_is_subn ? 1'b0 : 1'b1), c.sig };


      // debug
      // debug = ex[9];
      // debug = ~diff_base[9];
      debug = 0;

      // ==================== branch 1 ====================
      // 1.xxx(prod) +/- 0.xxx(c)
      begin : branch_1 //for multiply and align the prod and c, get the sum/diff, but not rounded or normalized

	 // count leading zeros
	 lzn = 0;		
	 for(int i=47;i>=25;i--)
	   if(prod[i]==1'b0)
	     lzn++;		// 0 ~ 23
	   else
	     break;

	 sft_bits = ~diff_base + 10'd1 - lzn; // the range of this variable ???

	 prod_sft = prod << lzn;
	 frac_c_sft = {frac_c, 24'b0} >> sft_bits;
	 if(sft_bits >= 48)
	   frac_c_sticky = |frac_c;
	 else if(sft_bits > 24)
	   frac_c_sticky = |(frac_c << (48 - sft_bits)); // 24 - (sft_bits-24)
	 else
	   frac_c_sticky = 0;

	 // expn = ex+2-lzn;	// ex+2-lzn =c.exp - diff_base - lzn
	 expn = c.exp - diff_base - lzn; // must be positive ??

	 is_true_sub = 1'b1;
	 
	 {carry, sum_diff_t, sum_diff_ext} = {prod_sft, 1'b0} + {is_true_sub, {frac_c_sft, frac_c_sticky}^{49{is_true_sub}}} + is_true_sub; // sum_diff is positive
	 
	 sum_diff = {carry, sum_diff_t[47:24], sum_diff_t[23:22], (|sum_diff_t[21:0] | sum_diff_ext)}; // sum, guard/round/sticky bit

	 // for lzd path (48 bits)
	 sum_diff_lzd = sum_diff_t;

	 // if(sft_bits == 1)	// need leading zeros detect
	 //   begin
	 //      $display("sum_diff[2:0] = %b", sum_diff[2:0]);
	 //      $display("frac_c_sft[23:0] = %b, sticky=%b", frac_c_sft[23:0], frac_c_sticky);
	 //      $display("sum_diff_t[23:0] = %b", sum_diff_t[23:0]);
	 //      $display("prod_sft = %b", sum_diff[2:0]);
	 //   end

 	 if(sum_diff[27]) begin	// the carry bit
	    round_bits = {sum_diff[3], |sum_diff[2:0]};
	    lsb_bit = sum_diff[4];
	 end
	 else if(sum_diff[26]) begin // the MSB
	    round_bits = {sum_diff[2], |sum_diff[1:0]};
	    lsb_bit = sum_diff[3];
	 end
	 else begin
	    round_bits = sum_diff[1:0];
	    lsb_bit = sum_diff[2];
	 end

	 // ==================================================
	 // now get round_bits, lsb_bit, sum_diff(28 bits), expn for the path which not need count leading zeros
	 // ==================================================

	 if(is_true_sub && frac_c_sft[47:46]==2'b01 && !sum_diff[26]) // or sft_bits<=1
	   begin
	      sel_lzd_path = 1;
	   end
	 else
	   begin
	      sel_lzd_path = 0;
	   end

	 // $display("sft_bits = %d", sft_bits);

	 if(~diff_base >= lzn)
	   debug = 1;	 
	 
      end // block: branch_1


      $display("sel_lzd_path = %d", sel_lzd_path);

      // ================================================================================
      // branch for round and normalize of the final sum/diff 
      // ================================================================================
      
      begin : branch_sum_lzd_path
	 bit [47:0] sum_diff_sft;
      	 bit 	    tmp_bit;
      	 bit [7:0]  lzd_expn;
	 bit [4:0]  lzn_t;
	 
      	 // count leading zeros (0 to 24)   --- 24 means the diff is exact 0
      	 lzn_t = 0;
      	 for(int i=46;i>=24;i--) // 0 ~ 23
	   if(sum_diff_lzd[i] == 1'b0)
      	     lzn_t++;
      	   else
      	     break;

	 // $display("sum_diff[2:0] = %b", sum_diff[2:0]);
	 
      	 if(expn > lzn_t+1)
      	   begin
	      sum_diff_sft = sum_diff_lzd << (lzn_t+1);	// --- first shift left for 1 bit ???
	      lzd_expn = expn - lzn_t - 1;
      	   end
      	 else
      	   begin
      	      sum_diff_sft = sum_diff_lzd << (expn-1);
      	      lzd_expn = 0;
      	   end

	 // may round cause the overflow ????
	 // {round_bit, sticky_bit} = {sum_diff_sft[23], |sum_diff_sft[22:0]};

	 
      	 // lzd_path_result.sign = frac_ext_borrow ^ exp_diff_borrow;
      	 lzd_path_result.sign = 0;
      	 lzd_path_result.exp = lzd_expn;
      	 lzd_path_result.sig = sum_diff_sft[46:24];

      	 lzd_path_excp = '0;

      end // block: branch_sum_lzd_path

      
	 
      begin
	 round_nearest_even = 0;
   	 unique case(fp_cfg.rounding_mode)
   	   RZ: round_inc_orig = 4'b0000;
   	   RN: begin
	      round_inc_orig = 4'b1000;
	      round_nearest_even = (round_bits==2'b10);
   	   end	
   	   RU: round_inc_orig = {4{~sign}}; // or result_sign
   	   RD: round_inc_orig = {4{sign}};
   	 endcase // unique case (fp_cfg.rounding_mode)
      end // always_comb

      round_inc_t = { (round_nearest_even ? lsb_bit : round_inc_orig[3]), round_inc_orig[2:0] };
   
      if(sum_diff[27]) begin
	 round_inc = round_inc_t;
      end
      else if(sum_diff[26]) begin
	 round_inc = round_inc_t>>1;
      end
      else begin
	 round_inc = round_inc_t>>2;
      end
      
      // assign overflow_val_sel = |round_inc_orig;
      // assign round_inc_t = { (round_nearest_even ? lsb_bit : round_inc_orig[3]), round_inc_orig[2:0] }; 

      sum_diff_rnd = sum_diff + round_inc;

      if(sum_diff_rnd[27])
	begin
	   sum_diff_rnd = sum_diff_rnd >> 1;
	   expn = expn + 1;
	end

      
   
   result.sign = 0;
   // result.exp = expn + sum_diff_rnd[27];
   result.exp = expn;
   // result.sig = sum_diff_ext[47:25];
   result.sig = sum_diff_rnd[25:3];


      if(sel_lzd_path)
	result = lzd_path_result;

      
      excp = '0;

   endfunction // float32_fma_new

   
//    function automatic void float32_fma_branch_1(
// 						output result,
// 						// for rnd path
// 						output bit [27:0] sum_diff,
// 						output bit [1:0] round_bits,
// 						output bit lsb_bit,
// 						// for lzd path
// 						output bit [47:0] sum_diff_lzd,
// 						output bit [7:0] expn,
						
// 						input bit [:] xx,
// 						input bit [9:0] diff_base,
// 						input bit [47:0] prod,
// 						input bit [23:0] frac_c,
// 						input bit [7:0] c_exp,
// 						input bit is_sub
// 						);

//       bit [4:0] 					  lzn;
//       bit [47:0] 					  frac_c_sft;
//       bit 						  frac_c_sticky;
      
//       bit [9:0] 					      sft_bits;	// the range is ok ??
//       bit [47:0] 					      prod_sft;
//       bit [47:0] 					      sum_diff_t;
//       bit 						      sum_diff_ext;
      
//       // bit [27:0] 					      sum_diff, sum_diff_rnd;
//       // bit 						      is_true_sub;

//       bit [7:0] 					      expn;
//       bit 						      carry;

//       // count leading zeros
//       lzn = 0;		
//       for(int i=47;i>=25;i--)
// 	if(prod[i]==1'b0)
// 	  lzn++;		// 0 ~ 23
// 	else
// 	  break;
      
//       sft_bits = ~diff_base + 10'd1 - lzn; // the range of this variable ???

//       prod_sft = prod << lzn;
//       frac_c_sft = {frac_c, 24'b0} >> sft_bits;
//       if(sft_bits >= 48)
// 	frac_c_sticky = |frac_c;
//       else if(sft_bits > 24)
// 	frac_c_sticky = |(frac_c << (48 - sft_bits)); // 24 - (sft_bits-24)
//       else
// 	frac_c_sticky = 0;

//       // expn = ex+2-lzn;	// ex+2-lzn =c.exp - diff_base - lzn
//       expn = c_exp - diff_base - lzn; // must be positive ??

//       {carry, sum_diff_t, sum_diff_ext} = {prod_sft, 1'b0} + {is_sub, {frac_c_sft, frac_c_sticky}^{49{is_sub}}} + is_sub; // sum_diff is positive

//       // sum/diff/lsb/round_bits for rnd_path
//       sum_diff = {carry, sum_diff_t[47:24], sum_diff_t[23:22], (|sum_diff_t[21:0] | sum_diff_ext)}; // sum, guard/round/sticky bit
      
//       if(sum_diff[27])		// the carry bit
// 	{lsb_bit, round_bits} = {sum_diff[4], sum_diff[3], |sum_diff[2:0]};
//       else if(sum_diff[26])	// the MSB
// 	{lsb_bit, round_bits} = {sum_diff[3], sum_diff[2], |sum_diff[1:0]};
//       else
// 	{lsb_bit, round_bits} = {sum_diff[2], sum_diff[1], sum_diff[0]};      
      
//       // sum/diff for lzd path
//       sum_diff_lzd = sum_diff_t; // carry/sum_diff_ext must be 0


//       if(is_true_sub && frac_c_sft[47:46]==2'b01 && !sum_diff_t[47]) // or sft_bits<=1
// 	begin
// 	   sel_lzd_path = 1;
// 	end
//       else
// 	begin
// 	   sel_lzd_path = 0;
// 	end


//       // return sum_diff if sft_bits is not 1
//       // return sum_diff_t 

      
//       if(sft_bits == 1)	// need leading zeros detect
// 	begin
// 	   $display("sum_diff[2:0] = %b", sum_diff[2:0]);
// 	   $display("frac_c_sft[23:0] = %b, sticky=%b", frac_c_sft[23:0], frac_c_sticky);	
// 	   $display("sum_diff_t[23:0] = %b", sum_diff_t[23:0]);
// 	   $display("prod_sft = %b", sum_diff[2:0]);
// 	end

// 	 // ==================================================
// 	 // now get round_bits, lsb_bit, sum_diff(28 bits), expn for the path which not need count leading zeros
// 	 // ==================================================


// 	 // $display("sft_bits = %d", sft_bits);


// 	 if(~diff_base >= lzn)
// 	   debug = 1;	 
	 
//    end // block: branch_1

// endfunction // float32_fma_branch_1

   
   
//    function automatic void float32_fma_new2(
// 				       output float32_t result,
// 				       output excp_t excp,
// 				       input fp_cfg_t fp_cfg,
// 				       input float32_t a, b, c,
// 				       // flags for debug
// 				       output bit debug
// 				       );
//       bit [9:0] 				  ex;
//       bit [47:0] 				  prod, prod_sft, prod_t;
//       bit [23:0] 				  frac_a, frac_b;
//       bit 					  sign;
//       bit [4:0] 				  lz_n;
//       // for rounding
//       bit 					  round, sticky;

//       bit 					  round_carry;
//       bit [1:0] 				  tmp_bits;
//       // after rounding
//       bit [23:0] 				  frac_prod_rnd;
//       bit 					  frac_rnd_carry;
//       // result
//       bit [7:0] 				  result_exp;
//       bit [22:0] 				  result_sig;


//       // variables for FMA
//       bit [9:0] 				  diff_base;
      
//       bit [23:0] 				  frac_c;

      
//       bit [2:0] 				  frac_c_grs; // g/r/s bits
//       bit 					  c_is_subn;
//       bit [47:0] 				  frac_c_ext; // if c shift right
      
//       bit [7:0] 				  expn; // intermediate exponent of the sum of FMA
//       bit [23:0] 				  sum;
//       bit [2:0] 				  sum_grs;
//       bit [23:0] 				  sum_rnd;
//       bit  					  sum_rnd_carry;
//       bit 					  carry;

//       bit [2:0] 				  prod_grs; // guard/round/sticky bits


//       // ==================== rounding (common for ??) ====================
//       bit [3:0] 				  round_inc_orig, round_inc_t, round_inc;
//       bit 					  lsb_bit;
//       bit [1:0] 				  round_bits;
//       bit 					  round_nearest_even;
//       // logic       overflow_val_sel;



//       // for branch 1
//       bit [4:0] 				  lzn;
	 
//       bit [47:0] 				  frac_c_sft;
//       bit 					  frac_c_sticky;
      
//       bit [9:0] 				  sft_bits;	// the range is ok ??
//       bit [47:0] 				  sum_diff_t;
//       bit 					  sum_diff_ext;

//       bit [27:0] 				  sum_diff, sum_diff_rnd;
      
//       bit 					  is_true_sub;

//       float32_t lzd_path_result;
//       excp_t lzd_path_excp;
      
//       bit 					  sel_lzd_path;

      
      

//       sign = a.sign ^ b.sign;	// neight a nor b is 0
      
//       frac_a = a.exp ? {1'b1, a.sig} : {a.sig, 1'b0};
//       frac_b = b.exp ? {1'b1, b.sig} : {b.sig, 1'b0};
//       prod = frac_a * frac_b;

//       ex = a.exp + b.exp - 127 - 1;

//       // for FMA
//       diff_base = c.exp - a.exp - b.exp + 126; //c.exp - ex - 2;
//       c_is_subn = c.exp == '0;
//       frac_c = { (c_is_subn ? 1'b0 : 1'b1), c.sig };


//       // debug
//       // debug = ex[9];
//       // debug = ~diff_base[9];
//       debug = 0;


//       // ==================== branch 1 ====================
//       // 1.xxx(prod) +/- 0.xxx(c)


//       $display("sel_lzd_path = %d", sel_lzd_path);

//       // ================================================================================
//       // branch for round and normalize of the final sum/diff 
//       // ================================================================================
      
//       begin : branch_sum_lzd_path
//       	 bit [22:0] sum_diff_sft;
//       	 bit 	    tmp_bit;
//       	 bit [7:0]  lzd_expn;
// 	 bit [4:0]  lzn_t;
	 
//       	 // count leading zeros (0 to 24)   --- 24 means the diff is exact 0
//       	 lzn_t = 0;
//       	 for(int i=25;i>=2;i--)
//       	   if(sum_diff[i] == 1'b0) // sum_diff is 28bits
//       	     lzn_t++;
//       	   else
//       	     break;


// 	 // $display("sum_diff[2:0] = %b", sum_diff[2:0]);

	 
//       	 if(expn > lzn_t+1)
//       	   begin
//       	      {sum_diff_sft, tmp_bit} = sum_diff[25:2] << (lzn_t+1); // 24 bits		--- first shift left for 1 bit ???
//       	      lzd_expn = expn - lzn_t - 1;
//       	   end
//       	 else
//       	   begin
//       	      {sum_diff_sft, tmp_bit} = sum_diff[25:2] << (expn-1);
//       	      lzd_expn = 0;
//       	   end
	 
//       	 // lzd_path_result.sign = frac_ext_borrow ^ exp_diff_borrow;
//       	 lzd_path_result.sign = 0;
//       	 lzd_path_result.exp = lzd_expn;
//       	 lzd_path_result.sig = sum_diff_sft;

//       	 lzd_path_excp = '0;
	   
//       end // block: branch_sum_lzd_path

	 
//       begin
// 	 round_nearest_even = 0;
//    	 unique case(fp_cfg.rounding_mode)
//    	   RZ: round_inc_orig = 4'b0000;
//    	   RN: begin
// 	      round_inc_orig = 4'b1000;
// 	      round_nearest_even = (round_bits==2'b10);
//    	   end	
//    	   RU: round_inc_orig = {4{~sign}}; // or result_sign
//    	   RD: round_inc_orig = {4{sign}};
//    	 endcase // unique case (fp_cfg.rounding_mode)
//       end // always_comb

//       round_inc_t = { (round_nearest_even ? lsb_bit : round_inc_orig[3]), round_inc_orig[2:0] };
   
//       if(sum_diff[27]) begin
// 	 round_inc = round_inc_t;
//       end
//       else if(sum_diff[26]) begin
// 	 round_inc = round_inc_t>>1;
//       end
//       else begin
// 	 round_inc = round_inc_t>>2;
//       end

      
//       // assign overflow_val_sel = |round_inc_orig;
//       // assign round_inc_t = { (round_nearest_even ? lsb_bit : round_inc_orig[3]), round_inc_orig[2:0] }; 

//       sum_diff_rnd = sum_diff + round_inc;

//       if(sum_diff_rnd[27])
// 	begin
// 	   sum_diff_rnd = sum_diff_rnd >> 1;
// 	   expn = expn + 1;
// 	end

      
   
//    result.sign = 0;
//    // result.exp = expn + sum_diff_rnd[27];
//    result.exp = expn;
//    // result.sig = sum_diff_ext[47:25];
//    result.sig = sum_diff_rnd[25:3];


//       if(sel_lzd_path)
// 	result = lzd_path_result;

      
//       excp = '0;

//    endfunction // float32_fma_new2
   
endpackage // float_arith
   
`endif //  `ifndef FLOAT_ARITH_PKG_SVH
   